function [X, Y, B] = modelxy(X1, Y1, X2, Y2)
    b = regress(Y1, [ones(numel(Y1), 1), X1], 0.05);

    d = Y2 - [ones(numel(Y2), 1), X2] * b;

    X = [X1; X2];
    Y = [Y1; Y2];
    s = numel(Y1) + 1;

    R = inf;
    B = [];
    D = [];
    for v = min(d):0.1:max(d)
        Y(s:end, :) = Y2 - v;

        [b, ~, r, ~, ~] = regress(Y, [ones(numel(Y), 1), X], 0.05);
        t = max(abs(r));

        if t < R
            R = t;
            B = b;
            D = v;
        end
    end

    Y(s:end, :) = Y2 - D;

    fprintf(1, '%f\n', R);
end
